///@file Triangle.cpp
///@brief Implementation of the class representing a triangle
///@author Arnaud Duval
///@version 0.0
///@date 2010/03/31


#include "../include/Triangle.h"
#include <cassert>
#include <cmath>

///@brief Default constructor
Triangle::Triangle()
{

}

///@brief Constructor with 3 given points
///@param DA point A
///@param DB point B
///@param DC point C
Triangle::Triangle(Point DA, Point DB, Point DC)
{
	A = DA;
	B = DB;
	C = DC;
}

///@brief Destructor
Triangle::~Triangle()
{

}

///@brief Return Area
///@big only wrk for 2D space (3D space is not yet implemented)
double Triangle::Area()
{
	assert(A.Dim() == 2);
	assert(B.Dim() == 2);
	assert(C.Dim() == 2);
	return 0.5 * fabs((B.Coord(1) - A.Coord(1))*(C.Coord(2) - A.Coord(2)) - (C.Coord(1)-A.Coord(1))*(B.Coord(2)-A.Coord(2)));
}

///@brief Return centroid
Point Triangle::Centroid()
{
	Point res(2);
	Point I(2);
	I.Coord(1) = (B.Coord(1)+C.Coord(1))/2.;
	I.Coord(2) = (B.Coord(2)+C.Coord(2))/2.;
	res.Coord(1) = A.Coord(1) + 2./3.*(I.Coord(1) - A.Coord(1));
	res.Coord(2) = A.Coord(2) + 2./3.*(I.Coord(2) - A.Coord(2));
	return res;
}

///@brief Return quadratic moment about axis x from Point O
double Triangle::Ixx_O()
{
	double A1, A2, A3, B1, B2, B3;
	double xA, xB, xC, yA, yB, yC;
	xA = A.Coord(1);
	xB = B.Coord(1);
	xC = C.Coord(1);
	yA = A.Coord(2);
	yB = B.Coord(2);
	yC = C.Coord(2);
	double res = 0.;


	if((fabs(xA - xB) > 1.E-6)&&(fabs(xB - xC) > 1.E-6)&&(fabs(xC - xA) > 1.E-6))
	{
		A1 = (yA-yB)/(xA-xB);
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += -pow(A1, 3.) * (pow(xA, 4.) - pow(xB, 4.)) / 12.;
		res += -pow(A1, 2.) * B1 * (pow(xA, 3.) - pow(xB, 3.)) / 3.;
		res += - A1 * pow(B1, 2.) * (pow(xA, 2.) - pow(xB, 2.)) / 2.;
		res += pow(A3, 3.) * (pow(xA, 4.) - pow(xB, 4.)) / 12.;
		res += pow(A3, 2.) * B3 * (pow(xA, 3.) - pow(xB, 3.)) / 3.;
		res += A3 * pow(B3, 2.) * (pow(xA, 2.) - pow(xB, 2.)) / 2.;
		res += -(pow(B1, 3.) - pow(B3, 3.)) * (xA - xB) / 3.;
		res += -pow(A2, 3.) * (pow(xB, 4.) - pow(xC, 4.)) / 12.;
		res += -pow(A2, 2.) * B2 * (pow(xB, 3.) - pow(xC, 3.)) / 3.;
		res += - A2 * pow(B2, 2.) * (pow(xB, 2.) - pow(xC, 2.)) / 2.;
		res += pow(A3, 3.) * (pow(xB, 4.) - pow(xC, 4.)) / 12.;
		res += pow(A3, 2.) * B3 * (pow(xB, 3.) - pow(xC, 3.)) / 3.;
		res += A3 * pow(B3, 2.) * (pow(xB, 2.) - pow(xC, 2.)) / 2.;
		res += -(pow(B2, 3.) - pow(B3, 3.)) * (xB - xC) / 3.;
	}
	else if(fabs(xA - xB) <= 1.E-6)
	{
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += -pow(A2, 3.) * (pow(xB, 4.) - pow(xC, 4.)) / 12.;
		res += -pow(A2, 2.) * B2 * (pow(xB, 3.) - pow(xC, 3.)) / 3.;
		res += - A2 * pow(B2, 2.) * (pow(xB, 2.) - pow(xC, 2.)) / 2.;
		res += pow(A3, 3.) * (pow(xB, 4.) - pow(xC, 4.)) / 12.;
		res += pow(A3, 2.) * B3 * (pow(xB, 3.) - pow(xC, 3.)) / 3.;
		res += A3 * pow(B3, 2.) * (pow(xB, 2.) - pow(xC, 2.)) / 2.;
		res += -(pow(B2, 3.) - pow(B3, 3.)) * (xB - xC) / 3.;
	}
	else if(fabs(xB - xC) <= 1.E-6)
	{
		A1 = (yA-yB)/(xA-xB);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += -pow(A1, 3.) * (pow(xA, 4.) - pow(xB, 4.)) / 12.;
		res += -pow(A1, 2.) * B1 * (pow(xA, 3.) - pow(xB, 3.)) / 3.;
		res += - A1 * pow(B1, 2.) * (pow(xA, 2.) - pow(xB, 2.)) / 2.;
		res += pow(A3, 3.) * (pow(xA, 4.) - pow(xB, 4.)) / 12.;
		res += pow(A3, 2.) * B3 * (pow(xA, 3.) - pow(xB, 3.)) / 3.;
		res += A3 * pow(B3, 2.) * (pow(xA, 2.) - pow(xB, 2.)) / 2.;
		res += -(pow(B1, 3.) - pow(B3, 3.)) * (xA - xB) / 3.;
	}
	else if(fabs(xC - xA) <= 1.E-6)
	{
		A2 = (yC-yB)/(xC-xB);
		A3 = (yA-yB)/(xA-xB);
		B2 = (xC*yB - xB*yC)/(xC-xB);
		B3 = (xA*yB - xB*yA)/(xA-xB);
		res += pow(A2, 3.) * (pow(xC, 4.) - pow(xB, 4.)) / 12.;
		res += pow(A2, 2.) * B2 * (pow(xC, 3.) - pow(xB, 3.)) / 3.;
		res +=  A2 * pow(B2, 2.) * (pow(xC, 2.) - pow(xB, 2.)) / 2.;
		res += -pow(A3, 3.) * (pow(xC, 4.) - pow(xB, 4.)) / 12.;
		res += -pow(A3, 2.) * B3 * (pow(xC, 3.) - pow(xB, 3.)) / 3.;
		res += -A3 * pow(B3, 2.) * (pow(xC, 2.) - pow(xB, 2.)) / 2.;
		res += (pow(B2, 3.) - pow(B3, 3.)) * (xC - xB) / 3.;
	}

	///@brief Points of the triangle must be provided in trigonometric sense

	return -res;
}

///@brief Return quadratic moment about axis y from Point O
double Triangle::Iyy_O()
{
	double A1, A2, A3, B1, B2, B3;
	double xA, xB, xC, yA, yB, yC;
	xA = A.Coord(1);
	xB = B.Coord(1);
	xC = C.Coord(1);
	yA = A.Coord(2);
	yB = B.Coord(2);
	yC = C.Coord(2);
	double res = 0.;


	if((fabs(xA - xB) > 1.E-6)&&(fabs(xB - xC) > 1.E-6)&&(fabs(xC - xA) > 1.E-6))
	{
		A1 = (yA-yB)/(xA-xB);
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += A1 * (pow(xB, 4.)/4. - pow(xA, 4.)/4.);
		res += A3 * (pow(xA, 4.)/4. - pow(xB, 4.)/4.);
		res += - ((B1 - B3)*(pow(xA, 3.) - pow(xB, 3.)))/3.;
		res += A2 * (pow(xC, 4.)/4. - pow(xB, 4.)/4.);
		res += A3 * (pow(xB, 4.)/4. - pow(xC, 4.)/4.);
		res += - ((B2 - B3)*(pow(xB, 3.) - pow(xC, 3.)))/3.;
	}
	else if(fabs(xA - xB) <= 1.E-6)
	{
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += A2 * (pow(xC, 4.)/4. - pow(xB, 4.)/4.);
		res += A3 * (pow(xB, 4.)/4. - pow(xC, 4.)/4.);
		res += - ((B2 - B3)*(pow(xB, 3.) - pow(xC, 3.)))/3.;
	}
	else if(fabs(xB - xC) <= 1.E-6)
	{
		A1 = (yA-yB)/(xA-xB);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += A1 * (pow(xB, 4.)/4. - pow(xA, 4.)/4.);
		res += A3 * (pow(xA, 4.)/4. - pow(xB, 4.)/4.);
		res += - ((B1 - B3)*(pow(xA, 3.) - pow(xB, 3.)))/3.;
	}
	else if(fabs(xC - xA) <= 1.E-6)
	{
		A1 = (yA-yB)/(xA-xB);
		A2 = (yB-yC)/(xB-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		res += A1 * (pow(xB, 4.)/4. - pow(xA, 4.)/4.);
		res += A2 * (pow(xA, 4.)/4. - pow(xB, 4.)/4.);
		res += - ((B1 - B2)*(pow(xA, 3.) - pow(xB, 3.)))/3.;
	}

	return -res;
}

///@brief Return product of inertia Ixy from point O
double Triangle::Ixy_O()
{
	double A1, A2, A3, B1, B2, B3;
	double xA, xB, xC, yA, yB, yC;
	xA = A.Coord(1);
	xB = B.Coord(1);
	xC = C.Coord(1);
	yA = A.Coord(2);
	yB = B.Coord(2);
	yC = C.Coord(2);
	double res = 0.;


	if((fabs(xA - xB) > 1.E-6)&&(fabs(xB - xC) > 1.E-6)&&(fabs(xC - xA) > 1.E-6))
	{
		A1 = (yA-yB)/(xA-xB);
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += (-1./24.)*(3. * pow(A1, 2.)*(pow(xA, 4.) - pow(xB, 4.)) + 8.* A1*B1*(pow(xA, 3.) - pow(xB, 3.)) - 3.*pow(A3, 2.)*(pow(xA, 4.) - pow(xB, 4.)) - 8.*A3*B3*(pow(xA, 3.) - pow(xB, 3.)) + 6.*(pow(B1, 2.) - pow(B3, 2.))*(pow(xA, 2.) - pow(xB, 2.)));
		res += (-1./24.)*(3. * pow(A2, 2.)*(pow(xB, 4.) - pow(xC, 4.)) + 8.* A2*B2*(pow(xB, 3.) - pow(xC, 3.)) - 3.*pow(A3, 2.)*(pow(xB, 4.) - pow(xC, 4.)) - 8.*A3*B3*(pow(xB, 3.) - pow(xC, 3.)) + 6.*(pow(B2, 2.) - pow(B3, 2.))*(pow(xB, 2.) - pow(xC, 2.)));
	}
	else if(fabs(xA - xB) <= 1.E-6)
	{
		A2 = (yB-yC)/(xB-xC);
		A3 = (yA-yC)/(xA-xC);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += (-1./24.)*(3. * pow(A2, 2.)*(pow(xB, 4.) - pow(xC, 4.)) + 8.* A2*B2*(pow(xB, 3.) - pow(xC, 3.)) - 3.*pow(A3, 2.)*(pow(xB, 4.) - pow(xC, 4.)) - 8.*A3*B3*(pow(xB, 3.) - pow(xC, 3.)) + 6.*(pow(B2, 2.) - pow(B3, 2.))*(pow(xB, 2.) - pow(xC, 2.)));
	}
	else if(fabs(xB - xC) <= 1.E-6)
	{
		A1 = (yA-yB)/(xA-xB);
		A3 = (yA-yC)/(xA-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B3 = (xA*yC - xC*yA)/(xA-xC);
		res += (-1./24.)*(3. * pow(A1, 2.)*(pow(xA, 4.) - pow(xB, 4.)) + 8.* A1*B1*(pow(xA, 3.) - pow(xB, 3.)) - 3.*pow(A3, 2.)*(pow(xA, 4.) - pow(xB, 4.)) - 8.*A3*B3*(pow(xA, 3.) - pow(xB, 3.)) + 6.*(pow(B1, 2.) - pow(B3, 2.))*(pow(xA, 2.) - pow(xB, 2.)));
	}
	else if(fabs(xC - xA) <= 1.E-6)
	{
		A1 = (yA-yB)/(xA-xB);
		A2 = (yB-yC)/(xB-xC);
		B1 = (xA*yB - xB*yA)/(xA-xB);
		B2 = (xB*yC - xC*yB)/(xB-xC);
		res += (-1./24.)*(3. * pow(A1, 2.)*(pow(xA, 4.) - pow(xB, 4.)) + 8.* A1*B1*(pow(xA, 3.) - pow(xB, 3.)) - 3.*pow(A2, 2.)*(pow(xA, 4.) - pow(xB, 4.)) - 8.*A2*B2*(pow(xA, 3.) - pow(xB, 3.)) + 6.*(pow(B1, 2.) - pow(B2, 2.))*(pow(xA, 2.) - pow(xB, 2.)));
	}

	return -res;
}